1 Effect of UPSTM-Based Decorrelation on Feature Discovery

1.0.1 Loading the libraries

library("FRESA.CAD")
library(readxl)
library(igraph)
library(umap)
library(tsne)
library(entropy)

op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)

1.1 Material and Methods

1.2 The Data

trainLabeled <- read.delim("~/GitHub/FCA/Data/trainSet.txt")
validLabeled <- read.delim("~/GitHub/FCA/Data/arcene_valid.txt")
wholeArceneSet <- rbind(trainLabeled,validLabeled)


wholeArceneSet$Labels <-  1*(wholeArceneSet$Labels > 0)
wholeArceneSet[,1:ncol(trainLabeled)] <- sapply(wholeArceneSet,as.double)

1.2.0.1 Standarize the names for the reporting

studyName <- "ARCENE"
dataframe <- wholeArceneSet
outcome <- "Labels"
thro <- 0.8
cexheat = 0.10

TopVariables <- 10

1.3 Generaring the report

1.3.1 Libraries

Some libraries

library(psych)
library(whitening)
library("vioplot")
library("rpart")

1.3.2 Data specs

pander::pander(c(rows=nrow(dataframe),col=ncol(dataframe)-1))
rows col
200 10000 
pander::pander(table(dataframe[,outcome]))
0 1
112 88 

varlist <- colnames(dataframe)
varlist <- varlist[varlist != outcome]

largeSet <- length(varlist) > 1500

1.3.3 Scaling the data

Scaling and removing near zero variance columns and highly co-linear(r>0.99999) columns


  ### Some global cleaning
  sdiszero <- apply(dataframe,2,sd) > 1.0e-16
  dataframe <- dataframe[,sdiszero]

  varlist <- colnames(dataframe)[colnames(dataframe) != outcome]
  tokeep <- c(as.character(correlated_Remove(dataframe,varlist,thr=0.99999)),outcome)
  dataframe <- dataframe[,tokeep]

  varlist <- colnames(dataframe)
  varlist <- varlist[varlist != outcome]
  
  iscontinous <- sapply(apply(dataframe,2,unique),length) >= 5 ## Only variables with enough samples



dataframeScaled <- FRESAScale(dataframe,method="OrderLogit")$scaledData

1.4 The heatmap of the data

numsub <- nrow(dataframe)
if (numsub > 1000) numsub <- 1000


if (!largeSet)
{

  hm <- heatMaps(data=dataframeScaled[1:numsub,],
                 Outcome=outcome,
                 Scale=TRUE,
                 hCluster = "row",
                 xlab="Feature",
                 ylab="Sample",
                 srtCol=45,
                 srtRow=45,
                 cexCol=cexheat,
                 cexRow=cexheat
                 )
  par(op)
}

1.4.0.1 Correlation Matrix of the Data

The heat map of the data


if (!largeSet)
{

  par(cex=0.6,cex.main=0.85,cex.axis=0.7)
  #cormat <- Rfast::cora(as.matrix(dataframe[,varlist]),large=TRUE)
  cormat <- cor(dataframe[,varlist],method="pearson")
  cormat[is.na(cormat)] <- 0
  gplots::heatmap.2(abs(cormat),
                    trace = "none",
  #                  scale = "row",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "Original Correlation",
                    cexRow = cexheat,
                    cexCol = cexheat,
                     srtCol=45,
                     srtRow=45,
                    key.title=NA,
                    key.xlab="|Pearson Correlation|",
                    xlab="Feature", ylab="Feature")
  diag(cormat) <- 0
  print(max(abs(cormat)))
}

1.5 The decorrelation


DEdataframe <- IDeA(dataframe,verbose=TRUE,thr=thro)
#> 
#>  Included: 7748 , Uni p: 0.01544774 , Uncorrelated Base: 920 , Outcome-Driven Size: 0 , Base Size: 920 
#> 
#> 
 1 <R=1.000,thr=0.900,N= 6808>, Top: 144( 307 ).[ 1 : 144 Fa= 144 : 0.900 ]( 144 , 2766 , 0 ),<|>Tot Used: 2910 , Added: 2766 , Zero Std: 0 , Max Cor: 1.000
#> 
 2 <R=1.000,thr=0.900,N= 6808>, Top: 429( 53 )....[ 1 : 429 Fa= 572 : 0.900 ]( 429 , 3267 , 144 ),<|>Tot Used: 4665 , Added: 3267 , Zero Std: 0 , Max Cor: 1.000
#> 
 3 <R=1.000,thr=0.900,N= 6808>, Top: 735( 85 ).......[ 1 : 735 Fa= 1305 : 0.900 ]( 734 , 3131 , 572 ),<|>Tot Used: 5607 , Added: 3131 , Zero Std: 0 , Max Cor: 1.000
#> 
 4 <R=1.000,thr=0.900,N= 6808>, Top: 705( 27 ).......=[ 2 : 705 Fa= 1996 : 0.910 ]( 704 , 2220 , 1305 ),<|>Tot Used: 6067 , Added: 2220 , Zero Std: 0 , Max Cor: 1.000
#> 
 5 <R=1.000,thr=0.900,N= 6808>, Top: 481( 27 )....[ 1 : 481 Fa= 2460 : 0.900 ]( 476 , 1361 , 1996 ),<|>Tot Used: 6447 , Added: 1361 , Zero Std: 0 , Max Cor: 0.999
#> 
 6 <R=0.999,thr=0.900,N= 6808>, Top: 302( 5 )...[ 1 : 302 Fa= 2746 : 0.900 ]( 299 , 742 , 2460 ),<|>Tot Used: 6684 , Added: 742 , Zero Std: 0 , Max Cor: 0.999
#> 
 7 <R=0.999,thr=0.900,N= 6808>, Top: 183( 1 ).=[ 2 : 183 Fa= 2907 : 0.918 ]( 178 , 317 , 2746 ),<|>Tot Used: 6735 , Added: 317 , Zero Std: 0 , Max Cor: 0.999
#> 
 8 <R=0.999,thr=0.900,N= 6808>, Top: 68( 4 )[ 1 : 68 Fa= 2959 : 0.910 ]( 66 , 99 , 2907 ),<|>Tot Used: 6744 , Added: 99 , Zero Std: 0 , Max Cor: 0.993
#> 
 9 <R=0.993,thr=0.900,N= 6808>, Top: 17( 2 )[ 1 : 17 Fa= 2968 : 0.900 ]( 14 , 23 , 2959 ),<|>Tot Used: 6747 , Added: 23 , Zero Std: 0 , Max Cor: 0.998
#> 
 10 <R=0.998,thr=0.900,N= 6808>, Top: 4( 1 )[ 1 : 4 Fa= 2971 : 0.900 ]( 4 , 5 , 2968 ),<|>Tot Used: 6747 , Added: 5 , Zero Std: 0 , Max Cor: 0.984
#> 
 11 <R=0.984,thr=0.900,N= 6808>, Top: 1( 1 )[ 1 : 1 Fa= 2972 : 0.900 ]( 1 , 1 , 2971 ),<|>Tot Used: 6747 , Added: 1 , Zero Std: 0 , Max Cor: 0.900
#> 
 12 <R=0.900,thr=0.800,N= 2704>, Top: 1025( 24 ).........[ 1 : 1025 Fa= 3314 : 0.884 ]( 973 , 1239 , 2972 ),<|>Tot Used: 6798 , Added: 1239 , Zero Std: 0 , Max Cor: 0.993
#> 
 13 <R=0.993,thr=0.900,N=  293>, Top: 141( 1 ).[ 1 : 141 Fa= 3346 : 0.900 ]( 138 , 144 , 3314 ),<|>Tot Used: 6798 , Added: 144 , Zero Std: 0 , Max Cor: 0.990
#> 
 14 <R=0.990,thr=0.900,N=  293>, Top: 35( 1 )[ 1 : 35 Fa= 3363 : 0.900 ]( 35 , 37 , 3346 ),<|>Tot Used: 6798 , Added: 37 , Zero Std: 0 , Max Cor: 0.913
#> 
 15 <R=0.913,thr=0.900,N=  293>, Top: 1( 1 )[ 1 : 1 Fa= 3364 : 0.900 ]( 1 , 1 , 3363 ),<|>Tot Used: 6798 , Added: 1 , Zero Std: 0 , Max Cor: 0.900
#> 
 16 <R=0.900,thr=0.800,N=  802>, Top: 231( 1 )..[ 1 : 231 Fa= 3402 : 0.800 ]( 223 , 308 , 3364 ),<|>Tot Used: 6803 , Added: 308 , Zero Std: 0 , Max Cor: 0.997
#> 
 17 <R=0.997,thr=0.900,N=   76>, Top: 38( 1 )[ 1 : 38 Fa= 3415 : 0.900 ]( 38 , 38 , 3402 ),<|>Tot Used: 6803 , Added: 38 , Zero Std: 0 , Max Cor: 0.972
#> 
 18 <R=0.972,thr=0.900,N=   76>, Top: 6( 1 )[ 1 : 6 Fa= 3419 : 0.900 ]( 6 , 7 , 3415 ),<|>Tot Used: 6803 , Added: 7 , Zero Std: 0 , Max Cor: 0.899
#> 
 19 <R=0.899,thr=0.800,N=  333>, Top: 72( 1 )[ 1 : 72 Fa= 3429 : 0.800 ]( 70 , 137 , 3419 ),<|>Tot Used: 6805 , Added: 137 , Zero Std: 0 , Max Cor: 0.978
#> 
 20 <R=0.978,thr=0.900,N=   24>, Top: 12( 1 )[ 1 : 12 Fa= 3431 : 0.900 ]( 12 , 12 , 3429 ),<|>Tot Used: 6805 , Added: 12 , Zero Std: 0 , Max Cor: 0.898
#> 
 21 <R=0.898,thr=0.800,N=  139>, Top: 45( 5 )[ 1 : 45 Fa= 3435 : 0.800 ]( 40 , 72 , 3431 ),<|>Tot Used: 6812 , Added: 72 , Zero Std: 0 , Max Cor: 0.984
#> 
 22 <R=0.984,thr=0.900,N=   18>, Top: 9( 1 )[ 1 : 9 Fa= 3436 : 0.900 ]( 9 , 9 , 3435 ),<|>Tot Used: 6812 , Added: 9 , Zero Std: 0 , Max Cor: 0.914
#> 
 23 <R=0.914,thr=0.900,N=   18>, Top: 1( 1 )[ 1 : 1 Fa= 3437 : 0.900 ]( 1 , 1 , 3436 ),<|>Tot Used: 6812 , Added: 1 , Zero Std: 0 , Max Cor: 0.883
#> 
 24 <R=0.883,thr=0.800,N=   37>, Top: 16( 1 )[ 1 : 16 Fa= 3438 : 0.800 ]( 16 , 20 , 3437 ),<|>Tot Used: 6813 , Added: 20 , Zero Std: 0 , Max Cor: 0.946
#> 
 25 <R=0.946,thr=0.900,N=    2>, Top: 1( 1 )[ 1 : 1 Fa= 3439 : 0.900 ]( 1 , 1 , 3438 ),<|>Tot Used: 6813 , Added: 1 , Zero Std: 0 , Max Cor: 0.893
#> 
 26 <R=0.893,thr=0.800,N=    7>, Top: 4( 1 )[ 1 : 4 Fa= 3440 : 0.800 ]( 3 , 3 , 3439 ),<|>Tot Used: 6813 , Added: 3 , Zero Std: 0 , Max Cor: 0.800
#> 
 27 <R=0.800,thr=0.800,N=    7>
#> 
 [ 27 ], 0.7999886 Decor Dimension: 6813 Nused: 6813 . Cor to Base: 3176 , ABase: 50 , Outcome Base: 0 
#> 
varlistc <- colnames(DEdataframe)[colnames(DEdataframe) != outcome]

pander::pander(sum(apply(dataframe[,varlist],2,var)))

63594442

pander::pander(sum(apply(DEdataframe[,varlistc],2,var)))

6756815

pander::pander(entropy(discretize(unlist(dataframe[,varlist]), 256)))

3.08

pander::pander(entropy(discretize(unlist(DEdataframe[,varlistc]), 256)))

1.63

1.5.1 The decorrelation matrix


if (!largeSet)
{

  par(cex=0.6,cex.main=0.85,cex.axis=0.7)
  
  UPSTM <- attr(DEdataframe,"UPSTM")
  
  gplots::heatmap.2(1.0*(abs(UPSTM)>0),
                    trace = "none",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "Decorrelation matrix",
                    cexRow = cexheat,
                    cexCol = cexheat,
                   srtCol=45,
                   srtRow=45,
                    key.title=NA,
                    key.xlab="|Beta|>0",
                    xlab="Output Feature", ylab="Input Feature")
  
  par(op)
}

1.6 The heatmap of the decorrelated data

if (!largeSet)
{

  hm <- heatMaps(data=DEdataframe[1:numsub,],
                 Outcome=outcome,
                 Scale=TRUE,
                 hCluster = "row",
                 cexRow = cexheat,
                 cexCol = cexheat,
                 srtCol=45,
                 srtRow=45,
                 xlab="Feature",
                 ylab="Sample")
  par(op)
}

1.7 The correlation matrix after decorrelation

if (!largeSet)
{

  cormat <- cor(DEdataframe[,varlistc],method="pearson")
  cormat[is.na(cormat)] <- 0
  
  gplots::heatmap.2(abs(cormat),
                    trace = "none",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "Correlation after IDeA",
                    cexRow = cexheat,
                    cexCol = cexheat,
                     srtCol=45,
                     srtRow=45,
                    key.title=NA,
                    key.xlab="|Pearson Correlation|",
                    xlab="Feature", ylab="Feature")
  
  par(op)
  diag(cormat) <- 0
  print(max(abs(cormat)))
}

1.8 U-MAP Visualization of features

1.8.1 The UMAP based on LASSO on Raw Data


if (nrow(dataframe) < 1000)
{
  classes <- unique(dataframe[1:numsub,outcome])
  raincolors <- rainbow(length(classes))
  names(raincolors) <- classes
  datasetframe.umap = umap(scale(dataframe[1:numsub,varlist]),n_components=2)
  plot(datasetframe.umap$layout,xlab="U1",ylab="U2",main="UMAP: Original",t='n')
  text(datasetframe.umap$layout,labels=dataframe[1:numsub,outcome],col=raincolors[dataframe[1:numsub,outcome]+1])
}

1.8.2 The decorralted UMAP

if (nrow(dataframe) < 1000)
{

  datasetframe.umap = umap(scale(DEdataframe[1:numsub,varlistc]),n_components=2)
  plot(datasetframe.umap$layout,xlab="U1",ylab="U2",main="UMAP: After IDeA",t='n')
  text(datasetframe.umap$layout,labels=DEdataframe[1:numsub,outcome],col=raincolors[DEdataframe[1:numsub,outcome]+1])
}

1.9 Univariate Analysis

1.9.1 Univariate



univarRAW <- uniRankVar(varlist,
               paste(outcome,"~1"),
               outcome,
               dataframe,
               rankingTest="AUC")

100 : V100 200 : V201 300 : V302 400 : V402 500 : V504
600 : V606 700 : V707 800 : V807 900 : V907 1000 : V1008
1100 : V1108 1200 : V1209 1300 : V1309 1400 : V1409 1500 : V1509
1600 : V1610 1700 : V1710 1800 : V1810 1900 : V1911 2000 : V2012
2100 : V2113 2200 : V2213 2300 : V2313 2400 : V2417 2500 : V2518
2600 : V2620 2700 : V2722 2800 : V2822 2900 : V2922 3000 : V3023
3100 : V3123 3200 : V3223 3300 : V3326 3400 : V3428 3500 : V3528
3600 : V3629 3700 : V3734 3800 : V3835 3900 : V3935 4000 : V4038
4100 : V4140 4200 : V4243 4300 : V4344 4400 : V4445 4500 : V4547
4600 : V4649 4700 : V4751 4800 : V4853 4900 : V4954 5000 : V5055
5100 : V5156 5200 : V5256 5300 : V5360 5400 : V5462 5500 : V5564
5600 : V5666 5700 : V5768 5800 : V5868 5900 : V5970 6000 : V6070
6100 : V6171 6200 : V6271 6300 : V6372 6400 : V6473 6500 : V6573
6600 : V6675 6700 : V6777 6800 : V6881 6900 : V6983 7000 : V7088
7100 : V7190 7200 : V7291 7300 : V7393 7400 : V7496 7500 : V7597
7600 : V7701 7700 : V7803 7800 : V7904 7900 : V8007 8000 : V8108
8100 : V8209 8200 : V8310 8300 : V8414 8400 : V8516 8500 : V8620
8600 : V8721 8700 : V8822 8800 : V8925 8900 : V9026 9000 : V9128
9100 : V9232 9200 : V9332 9300 : V9433 9400 : V9533 9500 : V9638
9600 : V9739 9700 : V9841 9800 : V9944




univarDe <- uniRankVar(varlistc,
               paste(outcome,"~1"),
               outcome,
               DEdataframe,
               rankingTest="AUC",
               )

100 : V100 200 : La_V201 300 : V302 400 : La_V402 500 : V504
600 : V606 700 : V707 800 : V807 900 : La_V907 1000 : La_V1008
1100 : La_V1108 1200 : V1209 1300 : La_V1309 1400 : La_V1409 1500 : La_V1509
1600 : La_V1610 1700 : La_V1710 1800 : La_V1810 1900 : La_V1911 2000 : La_V2012
2100 : La_V2113 2200 : La_V2213 2300 : V2313 2400 : La_V2417 2500 : La_V2518
2600 : V2620 2700 : La_V2722 2800 : V2822 2900 : V2922 3000 : La_V3023
3100 : La_V3123 3200 : La_V3223 3300 : V3326 3400 : V3428 3500 : La_V3528
3600 : La_V3629 3700 : La_V3734 3800 : La_V3835 3900 : La_V3935 4000 : La_V4038
4100 : La_V4140 4200 : La_V4243 4300 : V4344 4400 : V4445 4500 : La_V4547
4600 : La_V4649 4700 : La_V4751 4800 : V4853 4900 : La_V4954 5000 : La_V5055
5100 : V5156 5200 : La_V5256 5300 : V5360 5400 : La_V5462 5500 : La_V5564
5600 : La_V5666 5700 : La_V5768 5800 : La_V5868 5900 : La_V5970 6000 : La_V6070
6100 : La_V6171 6200 : La_V6271 6300 : La_V6372 6400 : La_V6473 6500 : V6573
6600 : La_V6675 6700 : La_V6777 6800 : V6881 6900 : La_V6983 7000 : La_V7088
7100 : La_V7190 7200 : V7291 7300 : La_V7393 7400 : V7496 7500 : La_V7597
7600 : La_V7701 7700 : La_V7803 7800 : La_V7904 7900 : V8007 8000 : La_V8108
8100 : La_V8209 8200 : La_V8310 8300 : La_V8414 8400 : La_V8516 8500 : V8620
8600 : La_V8721 8700 : La_V8822 8800 : La_V8925 8900 : La_V9026 9000 : V9128
9100 : La_V9232 9200 : La_V9332 9300 : La_V9433 9400 : V9533 9500 : La_V9638
9600 : La_V9739 9700 : La_V9841 9800 : V9944

1.9.2 Final Table


univariate_columns <- c("caseMean","caseStd","controlMean","controlStd","controlKSP","ROCAUC")

##top variables
topvar <- c(1:length(varlist)) <= TopVariables
tableRaw <- univarRAW$orderframe[topvar,univariate_columns]
pander::pander(tableRaw)
  caseMean caseStd controlMean controlStd controlKSP ROCAUC
V5005 314.7 72.9 239 83.6 0.18466 0.772
V4960 47.5 49.3 124 97.3 0.19534 0.751
V2309 43.1 45.3 113 89.0 0.18665 0.751
V8368 44.9 46.1 116 90.6 0.21091 0.751
V312 47.2 48.0 122 94.7 0.21678 0.750
V3365 46.3 46.9 119 92.5 0.22139 0.749
V9617 40.9 44.7 109 87.5 0.15591 0.749
V414 47.5 50.4 125 100.4 0.16265 0.749
V9092 33.9 63.1 124 132.6 0.00199 0.748
V1936 316.0 79.5 243 79.2 0.28495 0.748 


topLAvar <- univarDe$orderframe$Name[str_detect(univarDe$orderframe$Name,"La_")]
topLAvar <- unique(c(univarDe$orderframe$Name[topvar],topLAvar[1:as.integer(TopVariables/2)]))
finalTable <- univarDe$orderframe[topLAvar,univariate_columns]


pander::pander(finalTable)
  caseMean caseStd controlMean controlStd controlKSP ROCAUC
La_V7136 126.213 39.76 78.641 47.61 8.73e-01 0.795
La_V3417 7.380 10.46 -2.737 9.20 8.70e-05 0.779
La_V6838 11.780 6.29 2.208 11.57 3.11e-07 0.778
La_V2185 19.065 31.48 -6.343 25.57 1.54e-04 0.777
La_V571 -11.330 3.74 -7.865 3.90 1.51e-01 0.772
La_V2747 37.862 27.41 13.202 23.24 6.29e-03 0.766
La_V316 0.699 1.40 -0.609 1.64 5.75e-03 0.764
La_V4806 -3.900 6.34 13.100 26.09 4.26e-09 0.763
La_V3725 208.777 49.07 153.815 68.59 6.12e-01 0.752
La_V2243 9.055 11.91 -1.669 9.76 9.64e-01 0.748 

dc <- getLatentCoefficients(DEdataframe)
fscores <- attr(DEdataframe,"fscore")


pander::pander(c(mean=mean(sapply(dc,length)),total=length(dc),fraction=length(dc)/(ncol(dataframe)-1)))
mean total fraction
3.49 6540 0.664 

theCharformulas <- attr(dc,"LatentCharFormulas")


finalTable <- rbind(finalTable,tableRaw[topvar[!(topvar %in% topLAvar)],univariate_columns])


orgnamez <- rownames(finalTable)
orgnamez <- str_remove_all(orgnamez,"La_")
finalTable$RAWAUC <- univarRAW$orderframe[orgnamez,"ROCAUC"]
finalTable$DecorFormula <- theCharformulas[rownames(finalTable)]
finalTable$fscores <- fscores[rownames(finalTable)]

Final_Columns <- c("DecorFormula","caseMean","caseStd","controlMean","controlStd","controlKSP","ROCAUC","RAWAUC","fscores")

finalTable <- finalTable[order(-finalTable$ROCAUC),]
pander::pander(finalTable[,Final_Columns])
  DecorFormula caseMean caseStd controlMean controlStd controlKSP ROCAUC RAWAUC fscores
La_V7136 - (0.642)V6896 + V7136 126.213 39.76 78.641 47.61 8.73e-01 0.795 0.573 -1
La_V3417 - (0.946)V2076 - (0.127)V3374 + V3417 7.380 10.46 -2.737 9.20 8.70e-05 0.779 0.581 -2
La_V6838 + (0.013)V1474 - (1.016)V5947 + V6838 11.780 6.29 2.208 11.57 3.11e-07 0.778 0.670 0
La_V2185 - (0.723)V1380 + V2185 19.065 31.48 -6.343 25.57 1.54e-04 0.777 0.592 0
La_V571 - (0.829)V405 + V571 - (0.203)V7374 -11.330 3.74 -7.865 3.90 1.51e-01 0.772 0.598 -1
V5005 NA 314.739 72.85 239.304 83.62 1.85e-01 0.772 0.772 NA
La_V2747 + V2747 - (0.906)V3001 37.862 27.41 13.202 23.24 6.29e-03 0.766 0.624 3
La_V316 + V316 - (0.355)V1380 + (1.902)V6976 - (2.551)V8743 0.699 1.40 -0.609 1.64 5.75e-03 0.764 0.545 -3
La_V4806 + V4806 - (0.763)V7654 -3.900 6.34 13.100 26.09 4.26e-09 0.763 0.538 4
La_V3725 + V3725 - (0.542)V7868 208.777 49.07 153.815 68.59 6.12e-01 0.752 0.679 2
V4960 NA 47.511 49.25 123.696 97.33 1.95e-01 0.751 0.751 NA
V2309 NA 43.057 45.32 112.616 89.00 1.87e-01 0.751 0.751 NA
V8368 NA 44.886 46.06 116.036 90.62 2.11e-01 0.751 0.751 NA
V312 NA 47.182 48.01 121.634 94.73 2.17e-01 0.750 0.750 NA
V3365 NA 46.261 46.93 119.054 92.52 2.21e-01 0.749 0.749 NA
V9617 NA 40.886 44.75 108.848 87.49 1.56e-01 0.749 0.749 NA
V414 NA 47.466 50.45 125.348 100.36 1.63e-01 0.749 0.749 NA
V9092 NA 33.886 63.11 123.607 132.62 1.99e-03 0.748 0.748 NA
La_V2243 + V2243 - (0.869)V4001 - (0.198)V7347 9.055 11.91 -1.669 9.76 9.64e-01 0.748 0.611 -2
V1936 NA 315.955 79.48 243.062 79.21 2.85e-01 0.748 0.748 NA

1.10 Comparing IDeA vs PCA vs EFA

1.10.1 PCA

featuresnames <- colnames(dataframe)[colnames(dataframe) != outcome]
pc <- prcomp(dataframe[,iscontinous],center = TRUE,scale. = TRUE)   #principal components
predPCA <- predict(pc,dataframe[,iscontinous])
PCAdataframe <- as.data.frame(cbind(predPCA,dataframe[,!iscontinous]))
colnames(PCAdataframe) <- c(colnames(predPCA),colnames(dataframe)[!iscontinous]) 
#plot(PCAdataframe[,colnames(PCAdataframe)!=outcome],col=dataframe[,outcome],cex=0.65,cex.lab=0.5,cex.axis=0.75,cex.sub=0.5,cex.main=0.75)

#pander::pander(pc$rotation)


PCACor <- cor(PCAdataframe[,colnames(PCAdataframe) != outcome])


  gplots::heatmap.2(abs(PCACor),
                    trace = "none",
  #                  scale = "row",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "PCA Correlation",
                    cexRow = 0.5,
                    cexCol = 0.5,
                     srtCol=45,
                     srtRow= -45,
                    key.title=NA,
                    key.xlab="Pearson Correlation",
                    xlab="Feature", ylab="Feature")

1.10.2 EFA


EFAdataframe <- dataframeScaled

if (length(iscontinous) < 2000)
{
  topred <- min(length(iscontinous),nrow(dataframeScaled),ncol(predPCA)/2)
  if (topred < 2) topred <- 2
  
  uls <- fa(dataframeScaled[,iscontinous],nfactors=topred,rotate="varimax",warnings=FALSE)  # EFA analysis
  predEFA <- predict(uls,dataframeScaled[,iscontinous])
  EFAdataframe <- as.data.frame(cbind(predEFA,dataframeScaled[,!iscontinous]))
  colnames(EFAdataframe) <- c(colnames(predEFA),colnames(dataframeScaled)[!iscontinous]) 


  
  EFACor <- cor(EFAdataframe[,colnames(EFAdataframe) != outcome])
  
  
    gplots::heatmap.2(abs(EFACor),
                      trace = "none",
    #                  scale = "row",
                      mar = c(5,5),
                      col=rev(heat.colors(5)),
                      main = "EFA Correlation",
                      cexRow = 0.5,
                      cexCol = 0.5,
                       srtCol=45,
                       srtRow= -45,
                      key.title=NA,
                      key.xlab="Pearson Correlation",
                      xlab="Feature", ylab="Feature")
}

1.11 Effect on CAR modeling

par(op)
par(xpd = TRUE)
dataframe[,outcome] <- factor(dataframe[,outcome])
rawmodel <- rpart(paste(outcome,"~."),dataframe,control=rpart.control(maxdepth=3))
pr <- predict(rawmodel,dataframe,type = "class")

  ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
  if (length(unique(pr))>1)
  {
    plot(rawmodel,main="Raw",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
    text(rawmodel, use.n = TRUE,cex=0.75)
    ptab <- epiR::epi.tests(table(pr==0,dataframe[,outcome]==0))
  }


pander::pander(table(dataframe[,outcome],pr))
  0 1
0 109 3
1 13 75 
pander::pander(ptab$detail[c(5,3,4,6),])
  statistic est lower upper
5 diag.ac 0.920 0.873 0.954
3 se 0.852 0.761 0.919
4 sp 0.973 0.924 0.994
6 diag.or 209.615 57.738 761.000 

par(op)
par(xpd = TRUE)
DEdataframe[,outcome] <- factor(DEdataframe[,outcome])
IDeAmodel <- rpart(paste(outcome,"~."),DEdataframe,control=rpart.control(maxdepth=3))
pr <- predict(IDeAmodel,DEdataframe,type = "class")

  ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
  if (length(unique(pr))>1)
  {
    plot(IDeAmodel,main="IDeA",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
    text(IDeAmodel, use.n = TRUE,cex=0.75)
    ptab <- epiR::epi.tests(table(pr==0,DEdataframe[,outcome]==0))
  }

pander::pander(table(DEdataframe[,outcome],pr))
  0 1
0 108 4
1 11 77 
pander::pander(ptab$detail[c(5,3,4,6),])
  statistic est lower upper
5 diag.ac 0.925 0.879 0.957
3 se 0.875 0.787 0.936
4 sp 0.964 0.911 0.990
6 diag.or 189.000 58.011 615.766 

par(op)
par(xpd = TRUE)
PCAdataframe[,outcome] <- factor(PCAdataframe[,outcome])
PCAmodel <- rpart(paste(outcome,"~."),PCAdataframe,control=rpart.control(maxdepth=3))
pr <- predict(PCAmodel,PCAdataframe,type = "class")
ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
if (length(unique(pr))>1)
{
  plot(PCAmodel,main="PCA",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
  text(PCAmodel, use.n = TRUE,cex=0.75)
  ptab <- epiR::epi.tests(table(pr==0,PCAdataframe[,outcome]==0))
}

pander::pander(table(PCAdataframe[,outcome],pr))
  0 1
0 104 8
1 39 49 
pander::pander(ptab$detail[c(5,3,4,6),])
  statistic est lower upper
5 diag.ac 0.765 0.700 0.822
3 se 0.557 0.447 0.663
4 sp 0.929 0.864 0.969
6 diag.or 16.333 7.100 37.573 


par(op)

1.11.1 EFA


  EFAdataframe[,outcome] <- factor(EFAdataframe[,outcome])
  EFAmodel <- rpart(paste(outcome,"~."),EFAdataframe,control=rpart.control(maxdepth=3))
  pr <- predict(EFAmodel,EFAdataframe,type = "class")
  
  ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
  if (length(unique(pr))>1)
  {
    plot(EFAmodel,main="EFA",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
    text(EFAmodel, use.n = TRUE,cex=0.75)
    ptab <- epiR::epi.tests(table(pr==0,EFAdataframe[,outcome]==0))
  }


  pander::pander(table(EFAdataframe[,outcome],pr))
  0 1
0 109 3
1 13 75 
  pander::pander(ptab$detail[c(5,3,4,6),])
  statistic est lower upper
5 diag.ac 0.920 0.873 0.954
3 se 0.852 0.761 0.919
4 sp 0.973 0.924 0.994
6 diag.or 209.615 57.738 761.000 
  par(op)